data = readtable('../../data/main_VAR/US_full.xlsx');
date = data.year;

startdatenum = 1793;
enddatenum = 1946;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T_post = length(date(startdate:enddate));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option.figure = 0;

Tstart = startdate; 
Tend = enddate; 
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data((data.year >= startdatenum - 1), :);
data = data((data.year <= enddatenum), :);

taxrevgdp = data.tau(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.tau(1); % government tax revenue to GDP ratio
spendgdp = data.g(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.g(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.x(2:end); % real gdp growth (in logs)
inflation = data.pi(2:end); % growth in GDP price deflator (in logs)
ynom1 = data.x1yrCMT(2:end) / 100; % nominal yield on a 1-year Treasury bill, expressed per annum
ynom10 = data.x10yrCMT(2:end) / 100; % nominal yield on a 10-year Treasury bond, expressed per annum

ynom1 = log(ynom1 + 1);
ynom10 = log(ynom10 + 1);

pdm = data.pdm(2:end); % log price-dividend ratio on CRSP vw stock index; dividend is annual and seasonally adjusted
divgrm = data.divgrm(2:end); % nominal log dividend growth; dividend is annual and seasonally adjusted

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1;
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

% growth of log div to gdp ratio
divgrm = divgrm - inflation;

A0m = nanmean(pdm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tts(1) = taxrevgdp(1); 
gts(1) = spendgdp(1); 

for t = 2:T
    gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
    tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
end

disp('mean surplus')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Read in debt/gdp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gdebt = data.debt_gdp_sargenthall;

run ../../tools/run_var_estimate;

save(['MAT/res_annual_VAR_nodebt_US_before1946.mat'],'-regexp', '^(?!(option|Globaloption)$).');
